% MATLAB codes for Finite Element Analysis

%  Problem 1: 3 springs problem

% clear memory
clear ;

% ? define vars
P=10; % N
ki=[1;1;1;1];

% elementNodes: connections at elements
elementNodes=[1 2;
            2 3;
            2 4];

% numberElements: number of Elements
numberElements=size(elementNodes,1);

% numberNodes: number of nodes
numberNodes=4;

% for structure:
% displacements: displacement vector N x 1
% force : force vector N x 1
% stiffness: stiffness matrix N x N
displacements=zeros(numberNodes,1); %#ok<PREALL>
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes);

% applied load at node 2
force(2)=P;
% computation of the system stiffness matrix
for e=1:numberElements
    % elementDof: element degrees of freedom (Dof)
    elementDof=elementNodes(e,:) ;
    stiffness(elementDof,elementDof)=...
    stiffness(elementDof,elementDof)+ki(e)*[1 -1;-1 1];
end
% boundary conditions and solution
% prescribed dofs
prescribedDof=[1;3;4];
% free Dof : activeDof
% setdiff() 函数-->求两个数组的差集
% C = setdiff(A,B) 返回 A 中存在但 B 中不存在的数据,不包含重复项.C 是有序的.
activeDof=setdiff((1:numberNodes)',prescribedDof);

% solution
% !KU=F --> U=K\F
displacements=stiffness(activeDof,activeDof)\force(activeDof);

% positioning all displacements
displacements1=zeros(numberNodes,1);
displacements1(activeDof)=displacements;

% output displacements/reactions
outputDisplacementsReactions(displacements1,stiffness,numberNodes,prescribedDof)

